Brian E. J. Rose, University at Albany
You really should be looking at The Climate Laboratory book by Brian Rose, where all the same content (and more!) is kept up to date.
Here you are likely to find broken links and broken code.
This document uses the interactive Jupyter notebook format. The notes can be accessed in several different ways:
github at https://github.com/brian-rose/ClimateModeling_coursewareAlso here is a legacy version from 2015.
Many of these notes make use of the climlab package, available at https://github.com/brian-rose/climlab
I have run two sets of experiments with the CESM model:
Our main first task is to compute the two canonical measures of climate sensitivity for this model:
From the IPCC AR5 WG1 report, Chapter 9, page 817:
Equilibrium climate sensitivity (ECS) is the equilibrium change in global and annual mean surface air temperature after doubling the atmos- pheric concentration of CO2 relative to pre-industrial levels.
The transient climate response (TCR) is the change in global and annual mean surface temperature from an experiment in which the CO2 con- centration is increased by 1% yr–1, and calculated using the difference between the start of the experiment and a 20-year period centred on the time of CO2 doubling.
In [1]:
startingamount = 1.
amount = startingamount
for n in range(70):
amount *= 1.01
amount
Out[1]:
TCR is always smaller than ECS due to the transient effects of ocean heat uptake.
We are going to estimate the ECS of the fully coupled model by using the equilibrium response of the Slab Ocean .
In [2]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
In [3]:
casenames = {'cpl_control': 'cpl_1850_f19',
'cpl_CO2ramp': 'cpl_CO2ramp_f19',
'som_control': 'som_1850_f19',
'som_2xCO2': 'som_1850_2xCO2',
}
# The path to the THREDDS server, should work from anywhere
basepath = 'http://thredds.atmos.albany.edu:8080/thredds/dodsC/CESMA/'
# For better performance if you can access the roselab_rit filesystem (e.g. from JupyterHub)
#basepath = '/roselab_rit/cesm_archive/'
casepaths = {}
for name in casenames:
casepaths[name] = basepath + casenames[name] + '/concatenated/'
In [4]:
# make a dictionary of all the CAM atmosphere output
atm = {}
for name in casenames:
path = casepaths[name] + casenames[name] + '.cam.h0.nc'
print('Attempting to open the dataset ', path)
atm[name] = xr.open_dataset(path, decode_times=False)
In [5]:
days_per_year = 365
fig, ax = plt.subplots()
for name in ['cpl_control', 'cpl_CO2ramp']:
ax.plot(atm[name].time/days_per_year, atm[name].co2vmr*1E6, label=name)
ax.set_title('CO2 volume mixing ratio (CESM coupled simulations)')
ax.set_xlabel('Years')
ax.set_ylabel('pCO2 (ppm)')
ax.grid()
ax.legend();
Issues to think about:
The answer is closely related to the fact that the radiative forcing associated with CO2 increase is approximately logarithmic in CO2 amount. So a doubling of CO2 represents roughly the same radiative forcing regardless of the initial CO2 concentration.
In [6]:
# The surface air temperature, which we will use for our sensitivity metrics
atm['cpl_control'].TREFHT
Out[6]:
In [7]:
# The area weighting needed for global averaging
gw = atm['som_control'].gw
print(gw)
In [8]:
def global_mean(field, weight=gw):
'''Return the area-weighted global average of the input field'''
return (field*weight).mean(dim=('lat','lon'))/weight.mean(dim='lat')
In [9]:
# Loop through the four simulations and produce the global mean timeseries
TREFHT_global = {}
for name in casenames:
TREFHT_global[name] = global_mean(atm[name].TREFHT)
In [10]:
fig, axes = plt.subplots(2,1,figsize=(10,8))
for name in casenames:
if 'cpl' in name:
ax = axes[0]
ax.set_title('Fully coupled ocean')
else:
ax = axes[1]
ax.set_title('Slab ocean')
field = TREFHT_global[name]
field_running = field.rolling(time=12, center=True).mean()
line = ax.plot(field.time / days_per_year,
field,
label=name,
linewidth=0.75,
)
ax.plot(field_running.time / days_per_year,
field_running,
color=line[0].get_color(),
linewidth=2,
)
for ax in axes:
ax.legend();
ax.set_xlabel('Years')
ax.set_ylabel('Temperature (K)')
ax.grid();
ax.set_xlim(0,100)
fig.suptitle('Global mean surface air temperature in CESM simulations', fontsize=16);
Issues to think about here include:
In [11]:
# extract the last 10 years from the slab ocean control simulation
# and the last 20 years from the coupled control
nyears_slab = 10
nyears_cpl = 20
clim_slice_slab = slice(-(nyears_slab*12),None)
clim_slice_cpl = slice(-(nyears_cpl*12),None)
In [12]:
# extract the last 10 years from the slab ocean control simulation
T0_slab = TREFHT_global['som_control'].isel(time=clim_slice_slab).mean(dim='time')
T0_slab
Out[12]:
In [13]:
# and the last 20 years from the coupled control
T0_cpl = TREFHT_global['cpl_control'].isel(time=clim_slice_cpl).mean(dim='time')
T0_cpl
Out[13]:
In [14]:
# extract the last 10 years from the slab 2xCO2 simulation
T2x_slab = TREFHT_global['som_2xCO2'].isel(time=clim_slice_slab).mean(dim='time')
T2x_slab
Out[14]:
In [15]:
# extract the last 20 years from the coupled CO2 ramp simulation
T2x_cpl = TREFHT_global['cpl_CO2ramp'].isel(time=clim_slice_cpl).mean(dim='time')
T2x_cpl
Out[15]:
In [16]:
ECS = T2x_slab - T0_slab
TCR = T2x_cpl - T0_cpl
print('The Equilibrium Climate Sensitivity is {:.3} K.'.format(float(ECS)))
print('The Transient Climate Response is {:.3} K.'.format(float(TCR)))
Comparing against the multi-model mean of the ECS and TCR, our model is apparently slightly less sensitive than the CMIP5 mean.
Here is a helper function that takes a 2D lat/lon field and renders it as a nice contour map with accompanying zonal average line plot.
In [17]:
# The map projection capabilities come from the cartopy package. There are many possible projections
import cartopy.crs as ccrs
In [18]:
def make_map(field):
'''input field should be a 2D xarray.DataArray on a lat/lon grid.
Make a filled contour plot of the field, and a line plot of the zonal mean
'''
fig = plt.figure(figsize=(14,6))
nrows = 10; ncols = 3
mapax = plt.subplot2grid((nrows,ncols), (0,0), colspan=ncols-1, rowspan=nrows-1, projection=ccrs.Robinson())
barax = plt.subplot2grid((nrows,ncols), (nrows-1,0), colspan=ncols-1)
plotax = plt.subplot2grid((nrows,ncols), (0,ncols-1), rowspan=nrows-1)
cx = mapax.contourf(field.lon, field.lat, field, transform=ccrs.PlateCarree())
mapax.set_global(); mapax.coastlines();
plt.colorbar(cx, cax=barax, orientation='horizontal')
plotax.plot(field.mean(dim='lon'), field.lat)
plotax.set_ylabel('Latitude')
plotax.grid()
return fig, (mapax, plotax, barax), cx
In [19]:
# Plot a single time slice of surface air temperature just as example
fig, axes, cx = make_map(atm['cpl_control'].TREFHT.isel(time=0))
In [20]:
Tmap_cpl_2x = atm['cpl_CO2ramp'].TREFHT.isel(time=clim_slice_cpl).mean(dim='time')
Tmap_cpl_control = atm['cpl_control'].TREFHT.isel(time=clim_slice_cpl).mean(dim='time')
DeltaT_cpl = Tmap_cpl_2x - Tmap_cpl_control
Tmap_som_2x = atm['som_2xCO2'].TREFHT.isel(time=clim_slice_slab).mean(dim='time')
Tmap_som_control = atm['som_control'].TREFHT.isel(time=clim_slice_slab).mean(dim='time')
DeltaT_som = Tmap_som_2x - Tmap_som_control
In [21]:
fig, axes, cx = make_map(DeltaT_cpl)
fig.suptitle('Surface air temperature anomaly (coupled transient)', fontsize=16);
axes[1].set_xlim(0,7) # ensure the line plots have same axes
cx.set_clim([0, 8]) # ensure the contour maps have the same color intervals
fig, axes,cx = make_map(DeltaT_som)
fig.suptitle('Surface air temperature anomaly (equilibrium SOM)', fontsize=16);
axes[1].set_xlim(0,7)
cx.set_clim([0, 8])
Lots of intersting phenomena to think about here, including:
Continue to compare the transient and equilibrium responses in these CESM simulations.
Specifically, I would like you to examine the following:
Top-of-atmosphere energy budget:
The hydrological cycle:
Land - ocean warming contrast:
For all your results, please make an effort to point out any interesting or surprising results.
In [ ]:
In [22]:
%load_ext version_information
%version_information numpy, matplotlib, xarray, cartopy
Out[22]:
The author of this notebook is Brian E. J. Rose, University at Albany.
It was developed in support of ATM 623: Climate Modeling, a graduate-level course in the Department of Atmospheric and Envionmental Sciences
Development of these notes and the climlab software is partially supported by the National Science Foundation under award AGS-1455071 to Brian Rose. Any opinions, findings, conclusions or recommendations expressed here are mine and do not necessarily reflect the views of the National Science Foundation.
In [23]:
# # make a dictionary of all the CLM land model output
# land = {}
# for name in casenames:
# path = casepaths[name] + casenames[name] + '.clm2.h0.nc'
# print('Attempting to open the dataset ', path)
# land[name] = xr.open_dataset(path)
In [24]:
# # make a dictionary of all the sea ice model output
# ice = {}
# for name in casenames:
# path = casepaths[name] + casenames[name] + '.cice.h.nc'
# print('Attempting to open the dataset ', path)
# ice[name] = xr.open_dataset(path)
In [25]:
# # make a dictionary of all the river transport output
# rtm = {}
# for name in casenames:
# path = casepaths[name] + casenames[name] + '.rtm.h0.nc'
# print('Attempting to open the dataset ', path)
# rtm[name] = xr.open_dataset(path)
In [26]:
# ocn = {}
# for name in casenames:
# if 'cpl' in name:
# path = casepaths[name] + casenames[name] + '.pop.h.nc'
# print('Attempting to open the dataset ', path)
# ocn[name] = xr.open_dataset(path)
In [ ]: